#!/usr/bin/perl

require "tblas.pl";

package gamessio;

use feature 'switch';

BEGIN 
{
$bohr = 0.5291772108; # 1 a u = 0.5291772108 Angstrom
$to_kcal = 627.5; # 1 a u = 627.5 kkal/mole
$to_debye = 2.541922207165287; # 1 a u = 2.541922207165287 Debye
$to_buckingham = 1.345034181111014; # 1 a u = 1.345034181111014 Buckinghams

%atom_charges_names = 
(
"DUMMY"         =>   0, "HYDROGEN"      =>   1, "HELIUM"        =>   2,     
"LITHIUM"       =>   3, "BERYLLIUM"     =>   4, "BORON"         =>   5, "CARBON"        =>   6, "NITROGEN"      =>   7, "OXYGEN"        =>   8, "FLUORINE"      =>   9, "NEON"          =>  10, 
"SODIUM"        =>  11, "MAGNESIUM"     =>  12, "ALUMINIUM"     =>  13, "SILICON"       =>  14, "PHOSPHORUS"    =>  15, "SULFUR"        =>  16, "CHLORINE"      =>  17, "ARGON"         =>  18, 
"POTASSIUM"     =>  19, "CALCIUM"       =>  20, "SCANDIUM"      =>  21, "TITANIUM"      =>  22, "VANADIUM"      =>  23, "CHROMIUM"      =>  24, "MANGANESE"     =>  25, "IRON"          =>  26, "COBALT"        =>  27, "NICKEL"        =>  28, "COPPER"        =>  29, "ZINC"          =>  30, "GALLIUM"       =>  31, "GERMANIUM"     =>  32, "ARSENIC"       =>  33, "SELENIUM"      =>  34, "BROMINE"       =>  35, "KRYPTON"       =>  36, 
"RUBIDIUM"      =>  37, "STRONTIUM"     =>  38, "YTTRIUM"       =>  39, "ZIRCONIUM"     =>  40, "NIOBIUM"       =>  41, "MOLYBDENUM"    =>  42, "TECHNETIUM"    =>  43, "RUTHENIUM"     =>  44, "RHODIUM"       =>  45, "PALLADIUM"     =>  46, "SILVER"        =>  47, "CADMIUM"       =>  48, "INDIUM"        =>  49, "TIN"           =>  50, "ANTIMONY"      =>  51, "TELLURIUM"     =>  52, "IODINE"        =>  53, "XENON"         =>  54, 
"CAESIUM"       =>  55, "BARIUM"        =>  56, "LANTHANUM"     =>  57, "CERIUM"        =>  58, "PRASEODYMIUM"  =>  59, "NEODYMIUM"     =>  60, "PROMETHIUM"    =>  61, "SAMARIUM"      =>  62, "EUROPIUM"      =>  63, "GADOLINIUM"    =>  64, "TERBIUM"       =>  65, "DYSPROSIUM"    =>  66, "HOLMIUM"       =>  67, "ERBIUM"        =>  68, "THULIUM"       =>  69, "YTTERBIUM"     =>  70, "LUTETIUM"      =>  71, "HAFNIUM"       =>  72, "TANTALUM"      =>  73, "TUNGSTEN"      =>  74, "RHENIUM"       =>  75, "OSMIUM"        =>  76, "IRIDIUM"       =>  77, "PLATINUM"      =>  78, "GOLD"          =>  79, "MERCURY"       =>  80, "THALLIUM"      =>  81, "LEAD"          =>  82, "BISMUTH"       =>  83, "POLONIUM"      =>  84, "ASTATINE"      =>  85, "RADON"         =>  86, 
"FRANCIUM"      =>  87, "RADIUM"        =>  88, "ACTINIUM"      =>  89, "THORIUM"       =>  90, "PROTACTINIUM"  =>  91, "URANIUM"       =>  92, "NEPTUNIUM"     =>  93, "PLUTONIUM"     =>  94, "AMERICIUM"     =>  95, "CURIUM"        =>  96, "BERKELIUM"     =>  97, "CALIFORNIUM"   =>  98, "EINSTEINIUM"   =>  99, "FERMIUM"       => 100, "MENDELEVIUM"   => 101, "NOBELIUM"      => 102, "LAWRENCIUM"    => 103, "RUTHERFORDIUM" => 104, "DUBNIUM"       => 105, "SEABORGIUM"    => 106, "BOHRIUM"       => 107, "HASSIUM"       => 108, "MEITNERIUM"    => 109, "DARMSTADTIUM"  => 110, "ROENTGENIUM"   => 111, "UNUNBIUM"      => 112, "UNUNTRIUM"     => 113, "UNUNQUADIUM"   => 114, "UNUNPENTIUM"   => 115, "UNUNHEXIUM"    => 116
);                
                  
%atom_charges_symbols = 
(
"XX"  =>   0, "H"   =>   1, "HE"  =>   2,  
"LI"  =>   3, "BE"  =>   4, "B"   =>   5, "C"   =>   6, "N"   =>   7, "O"   =>   8, "F"   =>   9, "NE"  =>  10, 
"NA"  =>  11, "MG"  =>  12, "AL"  =>  13, "SI"  =>  14, "P"   =>  15, "S"   =>  16, "CL"  =>  17, "AR"  =>  18, 
"K"   =>  19, "CA"  =>  20, "SC"  =>  21, "TI"  =>  22, "V"   =>  23, "CR"  =>  24, "MN"  =>  25, "FE"  =>  26, "CO"  =>  27, "NI"  =>  28, "CU"  =>  29, "ZN"  =>  30, "GA"  =>  31, "GE"  =>  32, "AS"  =>  33, "SE"  =>  34, "BR"  =>  35, "KR"  =>  36, 
"RB"  =>  37, "SR"  =>  38, "Y"   =>  39, "ZR"  =>  40, "NB"  =>  41, "MO"  =>  42, "TC"  =>  43, "RU"  =>  44, "RH"  =>  45, "PD"  =>  46, "AG"  =>  47, "CD"  =>  48, "IN"  =>  49, "SN"  =>  50, "SB"  =>  51, "TE"  =>  52, "I"   =>  53, "XE"  =>  54, 
"CS"  =>  55, "BA"  =>  56, "LA"  =>  57, "CE"  =>  58, "PR"  =>  59, "ND"  =>  60, "PM"  =>  61, "SM"  =>  62, "EU"  =>  63, "GD"  =>  64, "TB"  =>  65, "DY"  =>  66, "HO"  =>  67, "ER"  =>  68, "TM"  =>  69, "YB"  =>  70, "LU"  =>  71, "HF"  =>  72, "TA"  =>  73, "W"   =>  74, "RE"  =>  75, "OS"  =>  76, "IR"  =>  77, "PT"  =>  78, "AU"  =>  79, "HG"  =>  80, "TL"  =>  81, "PB"  =>  82, "BI"  =>  83, "PO"  =>  84, "AT"  =>  85, "RN"  =>  86, 
"FR"  =>  87, "RA"  =>  88, "AC"  =>  89, "TH"  =>  90, "PA"  =>  91, "U"   =>  92, "NP"  =>  93, "PU"  =>  94, "AM"  =>  95, "CM"  =>  96, "BK"  =>  97, "CF"  =>  98, "ES"  =>  99, "FM"  => 100, "MD"  => 101, "NO"  => 102, "LR"  => 103, "RF"  => 104, "DB"  => 105, "SG"  => 106, "BH"  => 107, "HS"  => 108, "MT"  => 109, "DS"  => 110, "RG"  => 111, "UUB" => 112, "UUT" => 113, "UUQ" => 114, "UUP" => 115, "UUH" => 116
);

@atom_names = 
(
"Dummy"    , "Hydrogen" , "Helium"   ,     
"Lithium"  , "Beryllium", "Boron"    , "Carbon"   , "Nitrogen"    , "Oxygen"    , "Fluorine"  , "Neon"     , 
"Sodium"   , "Magnesium", "Aluminium", "Silicon"  , "Phosphorus"  , "Sulfur"    , "Chlorine"  , "Argon"    , 
"Potassium", "Calcium"  , "Scandium" , "Titanium" , "Vanadium"    , "Chromium"  , "Manganese" , "Iron"     , "Cobalt"   , "Nickel"    , "Copper"   , "Zinc"       , "Gallium"    , "Germanium", "Arsenic"    , "Selenium" , "Bromine"   , "Krypton"      , 
"Rubidium" , "Strontium", "Yttrium"  , "Zirconium", "Niobium"     , "Molybdenum", "Technetium", "Ruthenium", "Rhodium"  , "Palladium" , "Silver"   , "Cadmium"    , "Indium"     , "Tin"      , "Antimony"   , "Tellurium", "Iodine"    , "Xenon"        , 
"Caesium"  , "Barium"   , "Lanthanum", "Cerium"   , "Praseodymium", "Neodymium" , "Promethium", "Samarium" , "Europium" , "Gadolinium", "Terbium"  , "Dysprosium" , "Holmium"    , "Erbium"   , "Thulium"    , "Ytterbium", "Lutetium"  , "Hafnium"      , "Tantalum", "Tungsten"  , "Rhenium", "Osmium" , "Iridium"   , "Platinum"    , "Gold"       , "Mercury" , "Thallium" , "Lead"       , "Bismuth"    , "Polonium"  , "Astatine", "Radon", 
"Francium" , "Radium"   , "Actinium" , "Thorium"  , "Protactinium", "Uranium"   , "Neptunium" , "Plutonium", "Americium", "Curium"    , "Berkelium", "Californium", "Einsteinium", "Fermium"  , "Mendelevium", "Nobelium" , "Lawrencium", "Rutherfordium", "Dubnium" , "Seaborgium", "Bohrium", "Hassium", "Meitnerium", "Darmstadtium", "Roentgenium", "Ununbium", "Ununtrium", "Ununquadium", "Ununpentium", "Ununhexium"
);  
@atom_symbols = 
(
"Xx", "H" , "He",  
"Li", "Be", "B" , "C" , "N" , "O" , "F" , "Ne", 
"Na", "Mg", "Al", "Si", "P" , "S" , "Cl", "Ar", 
"K" , "Ca", "Sc", "Ti", "V" , "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", 
"Rb", "Sr", "Y" , "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I" , "Xe", 
"Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W" , "Re", "Os", "Ir", "Pt", "Au", "Hg" , "Tl" , "Pb" , "Bi" , "Po" , "At", "Rn", 
"Fr", "Ra", "Ac", "Th", "Pa", "U" , "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Uub", "Uut", "Uuq", "Uup", "Uuh"
);
@atom_vdw_radii =
(
0.00, 1.10, 
1.40, 1.81, 1.53, 1.92, 1.70, 1.55, 1.52, 1.47,
1.54, 2.27, 1.73, 1.84, 2.10, 1.80, 1.80, 1.75,
1.88, 2.75, 2.31, 2.30, 2.15, 2.05, 2.05, 2.05, 2.05, 2.00, 2.00, 2.00, 2.10, 1.87, 2.11, 1.85, 1.90, 1.83, 
2.02, 3.03, 2.49, 2.40, 2.30, 2.15, 2.10, 2.05, 2.05, 2.00, 2.05, 2.10, 2.20, 2.20, 1.93, 2.17, 2.06, 1.98, 
2.16, 3.43, 2.68, 2.50, 2.48, 2.47, 2.45, 2.43, 2.42, 2.40, 2.38, 2.37, 2.35, 2.33, 2.32, 2.30, 2.28, 2.27, 2.25, 2.20, 2.10, 2.05, 2.00, 2.00, 2.05, 2.10, 2.05, 1.96, 2.02, 2.07, 1.97, 2.02, 2.20,
3.48, 2.83, 2.00, 2.40, 2.00, 2.30, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 
);
@atom_masses =
(
0, 1.00794, 4.002602, 6.941, 9.012182, 10.811, 12.0107, 14.0067, 15.9994, 18.9984032, 20.1797, 22.98977, 24.3050, 26.981538, 28.0855, 30.973761, 32.065, 35.453, 39.948, 39.0983, 40.078, 44.95591, 47.867, 50.9415, 51.9961, 54.938049, 55.845, 58.9332, 58.6934, 63.546, 65.38, 69.723, 72.64, 74.92160, 78.96, 79.904, 83.798, 85.4678, 87.62, 88.90585, 91.224, 92.90638, 95.96, 98, 101.07, 102.90550, 106.42, 107.8682, 112.411, 114.818, 118.701, 121.760, 127.60, 126.90447, 131.293, 132.90545, 137.327, 138.9055, 140.116, 140.90765, 144.24, 145, 150.36, 151.964, 157.25, 158.92534, 162.500, 164.93032, 167.259, 168.93421, 173.054, 174.9668, 178.49, 180.9479, 183.84, 186.207, 190.23, 192.217, 195.078, 196.96655, 200.59, 204.3833, 207.2, 208.98040, 209, 210, 222, 223, 226, 227, 232.0381, 231.03588, 238.02891, 237, 244, 243, 247, 247, 251, 252, 257, 258, 259, 262, 267, 268, 271, 272, 270, 268, 281, 280, 285, 284, 289, 288, 293, 294, 294, 
);
  
};

sub check_normal_termination
# arguments: filename
# returns 1 when task terminated normally and 0 when it needs user's attention
{
  my $result = 0;
  open (my $out, '<', $_[0]) or die "out: $!";
  while (<$out>)
  {
    if (/NOT CONVERGED/)
    { 
      close $out or die "out: $!";
      return 0; 
    }
    if (/TERMINATED NORMALLY/) 
    { 
      close $out or die "out: $!";
      return 1; 
    }
  }       
}

sub read_eq_geometry_from_out
# arguments: filename
# returns: array: atom_name, charge, x, y, z
{
  my @result;
  open (my $out, '<', $_[0]) or die "out: $!";
  my $eq_geometry_header_found = 0;
  my $eq_geometry_found = 0;
  while (<$out>)
  {
    if ($eq_geometry_header_found == 1)
    {
      if ($eq_geometry_found == 1)
      {
        if (length($_) < 2)
        {
          $eq_geometry_found = 0;
          $eq_geometry_header_found = 0;
        } else {
          foreach ( split(/\s+/) ) { @result = (@result, $_) if (length($_) > 0); }          
        }
      } elsif ( /---------/ )
      {
        $eq_geometry_found = 1;
      }
    } elsif ( /EQUILIBRIUM GEOMETRY LOCATED/ ) {
      $eq_geometry_header_found = 1 
    }
  }
  close $out or die "out: $!";
  return @result;
}

sub read_geometry_from_inp
# arguments: filename
# returns: array: atom_name, charge, x, y, z
{
  my @result;
  open (my $out, '<', $_[0]) or die "out: $!";
  my $data_counter = 0;
  while (<$out>)
  {
    if ($data_counter != 0)
    {
      if (/\$END/) { $data_counter = 0; } elsif ( $data_counter > 2 )
      {
        my @tmp;
        foreach ( split(/\s+/) ) { @tmp = (@tmp, $_) if (length($_) > 0); }
        if ( $#tmp > 3 ) { push @result, @tmp; }  # to skip basis sets
      } else { ++$data_counter; }
    } elsif (/\$DATA/) { $data_counter = 1; }
  }
  close $out or die "out: $!";
  return @result;
}

sub print_geometry
# arguments: comment geometry
{
  my $result = " \$DATA\n".$_[0]."\n C1\n";
  my $n = (@_ - 1) / 5;
  for (my $i = 0; $i < $n; ++$i)
  {
    $result .= $_[1+$i*5]." ".$_[2+$i*5]." ".$_[3+$i*5]." ".$_[4+$i*5]." ".$_[5+$i*5]."\n";
  }
  return $result." \$END\n";
}

sub print_vdw_radii
# arguments: geometry
{
  my @geom = @_;
  my $result = " \$VDWRAD\n";
  for (my $i = 0; $i < @geom/5; ++$i)
  {
    $result .= " radius(".($i+1).")=".( ($geom[1+$i*5] > 0) ? $atom_vdw_radii[$geom[1+$i*5]] : 0 )."\n";
  }
  $result .= " \$END";
  return $result;
}

sub load_empty_ecp
# no args
{
  my %result;
  for (my $i = 0; $i < @atom_symbols; ++$i)
  {
    $result{uc($atom_symbols[$i])} = "  ".uc($atom_symbols[$i])." NONE\n";
  }
  return %result;
}

sub load_basis
# arguments filename, reference_to_wavefunctions_hash, reference_to_ecps_hash
# adds(and overrides) data to hashes
{
  my $wavefncref = $_[1];
  my $ecpref = $_[2];
  open (my $basislib, '<', $_[0]) or die "$!";
  my $current = "";
  my $wavefnc = 0;
  my $ecp = 0; 
  my $started = 0;
  while (<$basislib>)
  {
    if ($started == 1)
    {
     if ( /\$END/ ) { $started = 0; next; }
     my $line = $_;
     if ( /^!/ ) { next; }
     if ( /^([A-Za-z]{2,})/ )
     {
       if (length( $1) > 3 )
       {
         $current = uc $atom_symbols[ $atom_charges_names{$1} ];
       } else 
       {
         $current = $1;
       }
       $wavefnc = 1;
       $ecp = 0;
       $wavefncref->{$current} = "";
     } elsif ( /^  ([A-Za-z]{1,})-ECP/ )
     {
       $current = $1;
       $wavefnc = 0;
       $ecp = 1;
       $ecpref->{$current} = $line;
     } elsif ( length($current) > 0 ) 
     {
       if ($wavefnc == 1)
       {
         $wavefncref->{$current} .= $line;
       } elsif ($ecp == 1)
       {
         $ecpref->{$current} .= $line;
       }
     }
    } elsif ( /\$DATA/ || /\$ECP/) { $started = 1; } 
  }
}

sub print_geometry_with_basis
# arguments: comment, wavefunc_hash_ref, ecp_hash_ref, geometry   
{
  my $result = " \$DATA\n".$_[0]."\n C1\n";
  my $wavefncref = $_[1];
  my $ecpref = $_[2];
  my @geometry = @_[3..$#_];
  my $natoms = @geometry / 5;
  
  for (my $i = 0; $i < $natoms; ++$i )
  {
    $result .= $geometry[$i*5]." ".$geometry[1+$i*5]." ".$geometry[2+$i*5]." ".$geometry[3+$i*5]." ".$geometry[4+$i*5]."\n";
    $result .= $wavefncref->{uc($atom_symbols[$geometry[1+$i*5]])};
  }
  $result .= " \$END\n";
  my $ecp = " \$ECP\n";
  my $has_ecp = 0;
  for (my $i = 0; $i < $natoms; ++$i )
  {
    $ecp .= $ecpref->{uc($atom_symbols[$geometry[1+$i*5]])};
    if ( !( $ecpref->{uc($atom_symbols[$geometry[1+$i*5]])} =~ /NONE/ )) { $has_ecp = 1; }
  }
  $ecp .= " \$END\n"; 
  if ( $has_ecp == 1 )
  {
    $result .= $ecp;
  }
  return $result; 
}

sub generate_points_on_line
# arguments: center (3D vec) axis (3D vec) step npoints
{
  my @result;
  my @center = @_[0..2];
  my @axis = tblas::normalized_vector(@_[3..5]);
  my $step = $_[6];
  my $npoints = $_[7];
  my @point = tblas::sum_vectors(@center, tblas::scale( -$step * ($npoints + 1) / 2, @axis));
  my @vstep = tblas::scale( $step, @axis);
  for (my $i = 0; $i < $npoints; ++$i)
  {
    @point = tblas::sum_vectors(@point, @vstep);
    push @result, @point;
  }
  return @result;
}

sub print_points
# arguments: points
{
  my $n = @_ / 3;
  #print tblas::print_matrix(99, 3, @_);
  my $result = " \$POINTS\n ANGS ".$n."\n";
  for (my $i = 0; $i < $n; ++$i)
  {
    $result .= " ".$_[$i*3]." ".$_[$i*3+1]." ".$_[$i*3+2]."\n";
  }
  return $result." \$END\n";
}

sub print_grid
# arguments: step center xvec yvec xlen ylen
#
#             ylen
#  ............................
#  .......... xvec ............
#  ............^...............
#  ............|...............
#  .. yvec <-- 0 center ....... xlen
#  ............................
#  ............................
#  ............................
#  ............................
{
  my @center = @_[1..3];
  my @xvec = tblas::normalized_vector(@_[4..6]);
  my @yvec = tblas::normalized_vector(@_[7..9]);
  my @orign = tblas::sum_vectors( @center, tblas::sum_vectors( tblas::scale( -$_[10] / 2, @xvec), tblas::scale( -$_[11] / 2, @yvec) ) );
  my @vec1 = tblas::sum_vectors( @orign, tblas::scale($_[10], @xvec) );
  my @vec2 = tblas::sum_vectors( @orign, tblas::scale($_[11], @yvec) );
  my $result = " \$GRID UNITS=ANGS SIZE=".$_[0]."\n";
  $result .= " ORIGIN(1)=".$orign[0]." ORIGIN(2)=".$orign[1]."\n ORIGIN(3)=".$orign[2]."\n";
  $result .= " XVEC(1)=".$vec1[0]." XVEC(2)=".$vec1[1]." XVEC(3)=".$vec1[2]."\n";
  $result .= " YVEC(1)=".$vec2[0]." YVEC(2)=".$vec2[1]." YVEC(3)=".$vec2[2]."\n";
  $result .= " \$END\n";
  return $result;
}

sub read_potential_from_out
# arguments: filename
# returns array x y z val
# if potential was calculated in nucleus skips it
{
  my @result;
  open (my $out, '<', $_[0]) or die "out: $!";
  my $header_found = 0;
  my $data_counter = 0;
  my $skip = 0;
  while (<$out>)
  {
    if ( $header_found != 0 )
    {
      ++$data_counter;
      if ( $data_counter > 4 && length($_) > 2 )
      {
        if ( /WARNING/ ) 
        {
          $skip = 1;
        } elsif ( /END OF PROPERTY EVALUATION/ || /NET/)
        {
          $header_found = 0;
        } elsif ($skip != 1 && !/MERZ-KOLLMAN/ && !/VDW/ && !/NUMBER/) {
          push (@result, substr($_, 5, 9) * $bohr, substr($_, 15, 9) * $bohr, substr($_, 25, 9) * $bohr, substr($_, 61, 12));        
        } elsif ( !/MERZ-KOLLMAN/ && !/VDW/ && !/NUMBER/) { 
          $skip = 0;
          my @tmp;
          foreach ( split(/\s+/) ) { push (@tmp, $_) if (length($_) > 0); }
          push (@result, substr($_, 5, 9) * $bohr, substr($_, 15, 9) * $bohr, substr($_, 25, 9) * $bohr, 9999999999999);
        }
      }
    } elsif ( /ELECTROSTATIC POTENTIAL/ ) { $header_found = 1; @result = (); $data_counter = 0;  $skip = 0;  }
  }
  close $out or die "out: $!";
  return @result;
}

sub read_point_multipoles
#arguments: filename crd_array_ref charges_array_ref dipoles_array_ref quadrupoles_array_ref
#returns multipoles in a u coordinates in Angstroms
{
  my $filename = $_[0];
  my $crds = $_[1];
  my $chrgs = $_[2];
  my $dipoles = $_[3];
  my $quadrupoles = $_[4];
  
  open (my $out, '<', $filename) or die "out: $!";
  
  my $header_found = 0;
  my $counter = 0;
  while (<$out>)
  {
    my $line = $_;
    if ( /ELECTROSTATIC MOMENTS/ ) 
    {
      $header_found = 1;
      @$crds = ();
      @$chrgs = ();
      @$dipoles = ();
      @$quadrupoles = ();
    } elsif (/-----------------------/)
    {
      $header_found = 0;
      $counter = 0;
    } elsif ($header_found == 1) { 
      $counter++;
      if ($counter == 4)
      {
        my @tmp = split /\s+/, $line;
        push @$crds, $tmp[1]*$bohr, $tmp[2]*$bohr, $tmp[3]*$bohr;
        push @$chrgs, $tmp[4];
      } elsif ($counter == 6)
      {
        my @tmp = split /\s+/, $line;
        push @$dipoles, $tmp[1]/$to_debye, $tmp[2]/$to_debye, $tmp[3]/$to_debye;
      } elsif ($counter == 8)
      {
        my @tmp = split /\s+/, $line;
        push @$quadrupoles, $tmp[1]/$to_buckingham, $tmp[4]/$to_buckingham, $tmp[5]/$to_buckingham,
                            $tmp[4]/$to_buckingham, $tmp[2]/$to_buckingham, $tmp[6]/$to_buckingham,
                            $tmp[5]/$to_buckingham, $tmp[6]/$to_buckingham, $tmp[3]/$to_buckingham; 
      }
    }
  }
  
  close $out;
}

sub read_dma_multipoles
#arguments: filename crd_array_ref charges_array_ref dipoles_array_ref quadrupoles_array_ref
#returns multipoles in a u coordinates in Angstroms
{
  my $filename = $_[0];
  my $crds = $_[1];
  my $chrgs = $_[2];
  my $dipoles = $_[3];
  my $quadrupoles = $_[4];
  
  open (my $out, '<', $filename) or die "out: $!";
  
  my $header_found = 0;
  my $counter = 0;
  my @tmp;
  while (<$out>)
  {
    my $line = $_;
    if ( /DISTRIBUTED MULTIPOLAR ANALYSIS/ )
    {
      $header_found = 1; $counter = 0;
      @$crds = ();
      @$chrgs = ();
      @$dipoles = ();
      @$quadrupoles = ();
    } elsif ( /FIRST MOMENTS AT POINTS/ ) { $header_found = 2; $counter = 0; }
      elsif ( /SECOND MOMENTS AT POINTS/ ) { $header_found = 3; $counter = 0; }
      elsif ( /THIRD MOMENTS AT POINTS/ ) { $header_found = 0; $counter = 0; }
      elsif ( $header_found > 0 ) { $counter++; }
    if ( $counter >= 6 && $header_found == 1 )
    {
      if ( length > 5) 
      {
        push @tmp, $line; 
      } else {
        for (my $i = 0; $i < @tmp / 2; ++$i)
        {
          my @spl1;
          my @spl2;
          @spl1 = split /\s+/, $tmp[$i];
          @spl2 = split /\s+/, $tmp[$i+@tmp/2];
          push @$crds, $spl1[4]*$bohr, $spl1[5]*$bohr, $spl1[6]*$bohr;
          push @$chrgs, $spl1[3] + $spl2[3];
        }
        @tmp = ();
      }
    }
    if ( $counter >= 3 && $header_found == 2 )
    {
      if ( length > 5) 
      {
        @tmp = split /\s+/, $line;
        push @$dipoles, $tmp[2]/$to_debye, $tmp[3]/$to_debye, $tmp[4]/$to_debye;
        @tmp = ();
      } 
    }
    if ( $counter >= 3 && $header_found == 3 )
    {
      if ( length > 5) 
      {
        my @Q;
        @tmp = split /\s+/, $line;
        push @Q, $tmp[2]/$to_buckingham, $tmp[5]/$to_buckingham, $tmp[6]/$to_buckingham,
                 $tmp[5]/$to_buckingham, $tmp[3]/$to_buckingham, $tmp[7]/$to_buckingham,
                 $tmp[6]/$to_buckingham, $tmp[7]/$to_buckingham, $tmp[4]/$to_buckingham; 
        my $q=0.5*($Q[0]+$Q[4]+$Q[8]);
        push @$quadrupoles, 3/2*$Q[0]-$q, 3/2*$Q[1]   , 3/2*$Q[2]   ,
                            3/2*$Q[3]   , 3/2*$Q[4]-$q, 3/2*$Q[5]   ,
                            3/2*$Q[6]   , 3/2*$Q[7]   , 3/2*$Q[8]-$q;  
        @tmp = ();
      } 
    }

  }
}

sub calc_pot_charge
#arguments: distance charge in a u
#returns charge potential in a u
{
  return $_[1]/$_[0];
}

sub calc_pot_dipole
#arguments: axis distance dipole in a u
#returns dipole potential in a u
{
  return tblas::dot_prod(@_[0..2], @_[4..6])/$_[3]/$_[3];
}

sub calc_pot_quadrupole
#arguments: axis distance quadrupole in a u
#returns dipole potential in a u
{
  return tblas::dot_prod(@_[0..2], tblas::mult_matrix_vector_left(3, 3, @_[4..12], @_[0..2]))/$_[3]/$_[3]/$_[3];
}

sub calc_pot_multipoles
#arguments: point_crd_array crd_array charges_array dipoles_array quadrupoles_array
#returns multipoles potential in a u
{
  my @point = ($_[0]/$bohr, $_[1]/$bohr, $_[2]/$bohr);
  my @center = ($_[3]/$bohr, $_[4]/$bohr, $_[5]/$bohr);
  my $q = $_[6];
  my @dipole = @_[7..9];
  my @quadrupole = @_[10..18];
  
  my @axis = tblas::get_axis(@point, @center);
  my $dist = tblas::vector_length(tblas::sum_vectors(@point, tblas::scale(-1, @center)));
  
  my $result = $q/$dist + tblas::dot_prod(@axis, @dipole)/$dist/$dist + tblas::dot_prod(@axis, tblas::mult_matrix_vector_left(3, 3, @quadrupole, @axis))/$dist/$dist/$dist;

  return $result;
}

sub print_potential_esp
# arguments: q_tot, natoms, geometry, field
# returns a string formatted as esp
{
  my $result;
  my $natoms = $_[1];
  my $qtot = $_[0];
  my @geometry = @_[2..$natoms*5+1];
  my @pot = @_[$natoms*5+2..$#_];
  my $npoints = @pot / 4;
  
  $result = sprintf("%4d %5d %4d\n", $natoms, $npoints, $qtot);
  for (my $i = 0; $i < $natoms; ++$i)
  {
    $result .= sprintf("                 %15.7E %15.7E %15.7E\n", $geometry[2+$i*5]/$bohr, $geometry[3+$i*5]/$bohr, $geometry[4+$i*5]/$bohr);
  }
  for (my $i = 0; $i < $npoints; ++$i)
  {
    $result .= sprintf(" %15.7E %15.7E %15.7E %15.7E\n", $pot[3+$i*4], $pot[$i*4]/$bohr, $pot[$i*4+1]/$bohr, $pot[2+$i*4]/$bohr);
  }
  return $result;
}

sub check_line_on_field
# arguments: vec1 vec2 field (array of x y z val)
# returns 1 if vectors belongs to surface and 0 overwise
{
  my @vec1 = @_[0..2];
  my @vec2 = @_[3..5];
  my @pot = @_[6..$#_];
  my $step = tblas::distance(@pot[0..2], @pot[4..6]);
  my $eps = $step;
  my $npoints = @pot / 4;
  my $first = 0; my $second = 0;
  for (my $i = 0; $i < $npoints; ++$i)
  {
    if ( $first ==0 ) { if ( tblas::distance(@pot[$i*4..$i*4+2], @vec1) <= $eps ) { $first = 1; } } ;
    if ( $second ==0 ) { if ( tblas::distance(@pot[$i*4..$i*4+2], @vec2) <= $eps ) { $second = 1; } } ;
    if ( $first == 1 && $second == 1) { return 1; }
  }
  return 0;
}

sub interpolate_linear
{
# arguments vec1 val1 vec2 val2 vec3
# returns interpolated value in vec3
  my $val1 = $_[3];
  my $val2 = $_[7];
  my $dist1 = tblas::distance(@_[0..2], @_[8..10]);
  my $dist2 = tblas::distance(@_[4..6], @_[8..10]);
  return $val1 * $dist2 / ($dist1 + $dist2) + $val2 * $dist1 / ($dist1 + $dist2);
}

sub interpolate_bilinear
{
# arguments vec1 val1 vec2 val2 vec3 val3 vec4 val4 vec5 
# returns interpolated value in vec5
  my @vec1 = @_[0..2]; my $val1 = $_[3];
  my @vec2 = @_[4..6]; my $val2 = $_[7];
  my @vec3 = @_[8..10]; my $val3 = $_[11];
  my @vec4 = @_[12..14]; my $val4 = $_[15];
  my @vec5 = @_[16..18]; my $val5;
  
  if (tblas::compare_vectors(@vec5, @vec1) == 1 ) {return $val1;}
  if (tblas::compare_vectors(@vec5, @vec2) == 1 ) {return $val2;}
  if (tblas::compare_vectors(@vec5, @vec3) == 1 ) {return $val3;}
  if (tblas::compare_vectors(@vec5, @vec4) == 1 ) {return $val4;} 
  
  # sorting vertices
  my @v12 = tblas::get_axis(@vec1, @vec2);
  my @v13 = tblas::get_axis(@vec1, @vec3);
  my @v14 = tblas::get_axis(@vec1, @vec4);
  my @v23 = tblas::get_axis(@vec2, @vec3);
  my @v24 = tblas::get_axis(@vec2, @vec4);
  my @v34 = tblas::get_axis(@vec3, @vec4);
  if ( tblas::test_parallel(@v12, @v34) == 1 )
  {
    if (tblas::test_parallel( @v13, @v24 ) == 1 )
    {
      my @tmpvec = @vec1; @vec1 = @vec2; @vec2 = @tmpvec;
      my $tmpval = $val1; $val1 = $val2; $val2 = $tmpval;
    }
  } else {
    my @tmpvec = @vec2; @vec2 = @vec3; @vec3 = @tmpvec;
    my $tmpval = $val2; $val2 = $val3; $val3 = $tmpval;
  } 
  
  @v12 = tblas::sum_vectors(@vec2, tblas::scale(-1, @vec1));
  @v5 = tblas::sum_vectors(@vec5, tblas::scale(-1, @vec1));                      
  my $scale_factor = tblas::dot_prod(@v5, @v12)/tblas::dot_prod(@v12, @v12);
  
  my @p12 = tblas::sum_vectors(@vec1, tblas::scale($scale_factor, @v12));
  my @p34 = tblas::sum_vectors(@vec4, tblas::scale($scale_factor, @v12));
  
  return interpolate_linear( @p12, interpolate_linear(@vec1, $val1, @vec2, $val2, @p12), @p34, interpolate_linear(@vec3, $val3, @vec4, $val4, @p34), @vec5 );
}

sub get_field_x_y 
{
# returns number of points in field grid
# arguments field
  my $i = 1; 
  @field = @_;
  my $step_doubled = tblas::distance(@field[0..2], @field[4..6]) * 2;
  while ( tblas::distance(@field[$i*4-4..$i*4-2], @field[$i*4..$i*4+2]) < $step_doubled ) {++$i;}
  return ($i, @field / 4 / $i);
}

sub get_field_on_grid
{
# returns array x y z value
# arguments npoints points field_array
  my $npoints = $_[0];
  my @points = @_[1..$npoints*3];
  my @field = @_[1+$npoints*3..$#_];
  my $nfield = @field / 4;
  my @field_params = get_field_x_y(@field);
  my $fieldnx = $field_params[0];
  my $grid_step = tblas::distance( @field[0..2], @field[4..6] );
  
  my @xaxis = tblas::sum_vectors(@field[4..6], tblas::scale(-1, @field[0..2]));
  my @yaxis = tblas::sum_vectors(@field[4*$fieldnx..4*$fieldnx+2], tblas::scale(-1, @field[0..2]));
  my @orign = @field[0..2];
  
#  print "field: ".tblas::print_vector(@field_params)."\n";
   
  my @result;
  for (my $i = 0; $i < $npoints; ++$i)
  {
    @curpoint = tblas::sum_vectors(@points[$i*3..$i*3+2], tblas::scale(-1, @orign));
    my $n = (tblas::dot_prod(@curpoint, @xaxis)/tblas::dot_prod(@xaxis, @xaxis));
    my $m = (tblas::dot_prod(@curpoint, @yaxis)/tblas::dot_prod(@yaxis, @yaxis)); 
  #  print "n m ".$n." ".$m."\n";
    use POSIX;
    my $n1 = floor $n; my $n2 = ceil $n;
    my $m1 = floor $m; my $m2 = ceil $m;
  #  print tblas::print_vector($n1, $n2, $m1, $m2)."\n";
    push @result, @points[$i*3..$i*3+2], interpolate_bilinear( @field[($n1+$m1*$fieldnx)*4..($n1+$m1*$fieldnx)*4+3], @field[($n1+$m2*$fieldnx)*4..($n1+$m2*$fieldnx)*4+3], @field[($n2+$m2*$fieldnx)*4..($n2+$m2*$fieldnx)*4+3], @field[($n2+$m1*$fieldnx)*4..($n2+$m1*$fieldnx)*4+3], @points[$i*3..$i*3+2] );   
  }  
  return @result;
}

sub phantomize_atoms
{
#arguments natoms mask geometry
  my $natoms = $_[0];
  my @mask = @_[1..$natoms];
  my @geometry = @_[$natoms+1..$#_];
  my @result;
  for (my $i = 0; $i < $natoms; ++$i)
  {
    push @result, $geometry[$i*5], $geometry[$i*5 + 1] * $mask[$i], @geometry[$i*5+2..$i*5+4] ;
  }
  return @result;
}

return 1;

END {}